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Abstract 

We develop mathematical framework and computational tools for calculating frequency responses of linear 
time-invariant PDEs in which an independent spatial variable belongs to a compact interval. In conventional 
studies this computation is done numerically using spatial discretization of differential operators in the 
evolution equation. In this paper, we introduce an alternative method that avoids the need for finite- 
dimensional approximation of the underlying operators. This method recasts the frequency response operator 
as a two point boundary value problem and uses state-of-the-art automatic spectral collocation techniques for 
solving the resulting boundary value problems with accuracy comparable to machine precision. Our approach 
has two advantages over currently available schemes: first, it avoids numerical instabilities encountered in 
systems with differential operators of high order and, second, it alleviates difficulty in implementing boundary 
conditions. We provide examples from Newtonian and viscoelastic fluid dynamics to illustrate utility of the 
proposed method. 

Keywords: amplification of disturbances, automatic spectral collocation techniques, frequency responses, 
singular value decomposition, PDEs, spatio-temporal patterns, two point boundary value problems 



1. Introduction 

In many physical systems there is a need to examine the effects of exogenous disturbances on the variables 
of interest. The frequency response analysis represents an effective means for quantifying the system's 
performance in the presence of a stimulus, and it characterizes the steady-state response of a stable system 
to persistent harmonic forcing. At each temporal frequency, the frequency response of finite dimensional linear 
time-invariant systems with scalar input and output is a complex number that determines the magnitude and 
phase of the output relative to the input. In systems with many inputs and outputs (multi- variable systems), 
the frequency response is a complex matrix whose dimension is determined by the number of inputs and 
outputs. In systems with infinite dimensional input and output spaces, considered in this paper, the frequency 
response is an operator. It is well-known that the singular values of the frequency response matrix (in 
multi- variable systems) or the frequency response operator (in infinite dimensional systems) represent proper 
generalization of the magnitude characteristics for single-input single-output systems. At a specific frequency, 
the largest singular value determines the largest amplification from the input forcing to the desired output. 
Furthermore, the associated left and right principal singular functions identify the spatial distributions of 
the output (that exhibits this largest amplification) and the input (that has the strongest influence on the 
system's dynamics), respectively. 

In this paper, we study the frequency responses of linear time-invariant partial differential equations 
(PDEs) in which an independent spatial variable belongs to a compact interval. We are interested in com- 
puting the largest singular value of the frequency response operator and its corresponding singular functions. 
Computation of frequency responses for PDEs is typically done numerically using finite-dimensional approx- 
imations of the operators in the evolution equation. Pseudo-spectral methods represent a powerful tool for 
discretization of spatial differential operators, and they possess superior numerical accuracy compared to 
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approximation schemes based on finite differences [TJ [5J [31 H] . In spite of their advantages, pseudo-spectral 
methods may produce unreliable results and even fail to converge upon grid refinement when dealing with 
systems that contain differential operators of high order; this lack of convergence is attributed to the loss of 
accuracy arising from ill-conditioning of the discretized differentiation matrices [5]. Furthermore, implemen- 
tation of general boundary conditions may be challenging. 

To alleviate these difficulties, we introduce a method that avoids the need for finite dimensional approxi- 
mations of differential operators in the evolution equation. This is accomplished by recasting the frequency 
response operator as a two point boundary value problem (TPBVP) that is given by either an input-output 
differential equation of a high order or by an equivalent system of first order differential equations (i.e., 
spatial state-space representation). Furthermore, we present a procedure for converting these differential 
representations into the corresponding systems of integral equations. This transformation facilitates the use 
of recently developed computing environment, Chebfun 6J, that is capable of solving boundary value prob- 
lems and eigenvalue problems with superior accuracy. Our mathematical framework in conjunction with 
Chebfun's state-of-the-art numerical algorithms has two main advantages over standard methods: first, it 
alleviates numerical ill-conditioning encountered in systems with differential operators of high order; and 
second, it enables easy implementation of a wide range of boundary conditions. 

Chebfun is a collection of powerful algorithms for numerical computations that involve continuous and 
piecewise-continuous functions. Instead of working in a finite dimensional setting, Chebfun allows users to 
symbolically represent functions and operators in their infinite dimensional form with simple and compact 
Matlab syntaxes. This provides an elegant high-level language for solving linear and nonlinear boundary 
value and eigenvalue problems with few lines of code. Internally, functions are computed numerically us- 
ing automatic Chebyshev polynomial interpolation techniques, and the operators are approximated using 
automatic spectral collocation methods. Finite dimensional approximations of functions and operators are 
automatically refined in order to obtain accurate and convergent representations. Furthermore, once the 
boundary conditions are specified Chebfun makes sure that they are automatically satisfied internally when 
solving differential or integral equations. 

The proposed method has many potential applications in numerical analysis, physics, and engineering, 
especially in systems with generators that do not commute with their adjoints [Tj. In these systems, standard 
modal analysis may fail to capture amplification of exogenous disturbances, low stability margins, and large 
transient responses. In contrast, singular value decomposition of the frequency response operator represents 
an effective tool for identifying these non-modal aspects of the system's dynamics. In particular, wall-bounded 
shear flows of both Newtonian and viscoelastic fluids have non-normal dynamical generators of high spatial 
order and the ability to accurately compute frequency responses for these systems is of paramount importance; 
additional examples of systems with non-normal generators, for which the tools developed in this papers are 
particularly well-suited, can be found in the outstanding book by Trefethen and Embree [7J and the references 
therein. The utility of non-modal analysis in understanding the dynamics of infinitesimal fluctuations around 
laminar flow conditions has been well-documented; see El [lOj QT] for Newtonian fluids and (TSJ [T3J HU [15] 
for viscoelastic fluids. In viscoelastic fluids with large polymer relaxation times, analysis is additionally 
complicated by the fact that pseudo-spectral methods exhibit spurious numerical instabilities [TU [TTJ . We 
use examples from fluid mechanics to demonstrate the ease of incorporating boundary conditions and superior 
accuracy of our method compared to conventional finite dimensional approximation schemes. 

Our presentation is organized as follows. In Section [2j we formulate the problem and discuss the notion 
of a frequency response for PDEs. In Section [3j we present the method for converting the frequency response 
operator into a TPBVP that can be posed as an input-output differential equation or as a spatial state-space 
representation. In Section [4j we show how to transform a family of differential equations into equivalent 
integral equations and describe the use of Chebfun's indefinite integration operator for determining the eigen- 
values and corresponding eigenfunctions of the resulting integral equations. In Section [5] we demonstrate 
the utility of our developments by providing two examples from Newtonian and viscoelastic fluid dynamics. 
We conclude with a brief summary of the paper in Section |6j and relegate the mathematical developments 
to the appendices. 
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2. Motivating examples and problem formulation 



In this section, we provide two examples that are used to motivate our developments and to illustrate the 
classes of systems that we consider. These examples are used throughout the paper to explain the problem 
setup and utility of the proposed method. We then describe the class of PDEs that we study and briefly 
review the notion of a frequency response operator. 

2.1. Motivating examples 

We next present two physical examples: the one-dimensional (ID) diffusion equation, and the system of 
PDEs that governs the dynamics of the flow fluctuations in an inertialess channel flow of viscoclastic fluids. 
The ID diffusion equation has simple dynamics and it is used to illustrate mathematical framework developed 
in the paper. The example from viscoelastic fluid mechanics is used to demonstrate utility of our approach 
on a system that is known to produce spurious numerical instabilities. We show how numerical difficulties 
encountered in the computation of the frequency responses can be overcome using the developed framework 
in conjunction with state-of-the-art automatic spectral collocation techniques. 

2.1.1. One- dimensional diffusion equation 

Let a one-dimensional diffusion equation with homogenous Dirichlct boundary conditions and zero initial 
conditions be subject to spatially and temporally distributed forcing d(y,t), 

<t>t{y,t) = <t>w(y,t) + d(y,t), 

0(±M) = 0, (1) 
<b(y, 0) = 0, ye[-l, 1]. 

Throughout the paper, the spatially independent variable is denoted by y, the time is denoted by t, and 
the subscripts denote differentiation with respect to time/space. Considering (j> as the field of interest, the 
frequency response operator for this system (from input d to output <fi) is given by 

7» = (iul - D {2) y\ (2) 

where is the second derivative operator with homogenous Dirichlet boundary conditions, / is the identity 
operator, u> is the temporal frequency, and i is the imaginary unit. 

It is well known that the second derivative operator with Dirichlet boundary conditions is self-adjoint with 
a complete set of orthonormal eigenfunctions, v n (y) — sin {(niv/2)(y + 1)), n = {1, 2, . . .}. This information 
can be used to diagonalize operator in T{lo) which facilitates straightforward determination of the 

frequency response. For systems with spatially varying coefficients and non-normal generators the frequency 
response analysis is typically done numerically using finite dimensional approximations of the differential 
operators. For example, the pseudospectral method [18] with N collocation points can be used to transform 
the frequency response operator ^ of system (JlJ into an N x N matrix. However, for systems with differential 
operators of high order, spectral differentiation matrices may be poorly conditioned and implementation of 
boundary conditions may be challenging. In this paper, numerical approximation of differential operators in 
the evolution equation is avoided by first recasting the system into corresponding two point boundary value 
problems and then using state-of-the-art techniques for solving the resulting boundary value problems with 
accuracy comparable to machine precision. 

Alternatively, by applying the temporal Fourier transform on system ([T]) we obtain the following input- 
output differential equation 

&"(y,u) - iui(f>(y,uj) = -d(y,w), (3a) 
4>(±l,u) = 0, (3b) 



where d and <f> are the Fourier transformed input and output fields, and 4>' = &4>/&y. At each oj, (3a 



second-order ordinary differential equation (in y) subject to the boundary conditions (3b). Equivalently, by 
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Figure 1: We consider the dynamics of flow fluctuations in the (x,j/)-plane of the channel. 
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We will utilize structures of the TPBVPs ^ and Q in conjunction with recently developed automatic 
spectral collocation techniques to study the frequency response across u>. 

2.1.2. Inertialess channel flow of viscoelastic fluids 

We next consider a system that describes the dynamics of two-dimensional velocity and polymer stress 
fluctuations in an inertialess channel flow of viscoelastic fluids; see figure [I] for geometry. The dynamics of 
infinitesimal fluctuations around the mean flow (v, r) are given by 
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v = r u v ] T , p 1 and t arc the velocity, pressure, and stress fluctuations; u and v are velocities in x and y 



directions; V is the gradient; and V = V • V is the Laplacian. System (|5j) is driven by spatially distributed 
and temporally varying body force fluctuations d = [ di efe ] with d\ and c?2 representing the forcing in 
x and y. The non-dimensional parameters in ^ are the ratio of the solvent to the total viscosity /? S (0, 1), 
and the ratio of the fluid relaxation time to the characteristic flow time We (the Weissenberg number). 



Static- in-time momentum ( 5a ) and continuity ([5b]) equations describe the motion of incompressible fluids 
in the Stokes flow, i.e., at zero Reynolds number. The constitutive equation ( |5c[ ) captures the influence of 
the velocity gradients on the dynamics of stress fluctuations in dilute polymer solutions [T!5] • For background 
material on the use of frequency response analysis in understanding the dynamics of viscoelastic fluids, we 
refer the reader to 021 dU [H] and the references therein. 

By expressing the velocity fluctuations in terms of the stream function ip, 



d y ip, v 



-d x ip, 



4 



pressure can be removed form the equations (p)]). Furthermore, by applying the Fourier transform in x and 
t on (5c) and by substituting the resulting expression for stresses into the equation for ip, we arrive at the 
following ordinary differential equation (ODE) in y for the stream function, 



(D^ + a 3 (y)DW + a 2 (y)DW + a x {y)D^ + a (y)) 4>(y) = 
(b 1 (y)D( 1 )+b (y))d( 2/ ), 













ifca; 



T(w) 

to) 

= 4>(y = ±l) = ^'(y = ±l), 
where — d k /dy k , k x is the horizontal wavenumber, and 

DW o 

D (1) = 

DW 



(0) 



The coefficients bj (?/)} in ([6]) are reported in Appendix D The system of equations ([6]) is parameterized 

by u, fcj,, /3, and We. For notational convenience, we have suppressed the dependence of ip, d, u, and v on 
these four parameters. 

In Section [5] we show that spatial discretization of the operators in ([5| using the pseudo-spectral 
method [18j can produce erroneous frequency responses. In contrast, transformation of the system into 
a TPBVP followed by the use of Chebfun environment [6] yields reliable results. 

2.2. Problem formulation 

We now formulate the problem for PDEs with the evolution equation 



£ 0t(j/,*) = F<t>{y,t) - 
(p(y,t) = H4>{y,t), 



Qd(y,t), 



(7a) 
(7b) 



where t € [0,oo) and y £ [a,b] denote the temporal and spatial variables. The spatially distributed and 
temporally varying state, input, and output fields are represented by </>, d, and ip, respectively. At each 
t, d(-,i) and <p(-,t) denote the square-integrable vector-valued functions, and {£, Q, H} are matrices 
of differential operators with, in general, spatially varying coefficients. For example, the ijth entry of the 
operator J- can be expressed as 



k = 



where each is a function that is at least k times continuously differentiable on the interval [a, b] |20j 
]j( k ) — d k /dy k , and ray is the order of the highest derivative in .Fy. 

The application of the temporal Fourier transform yields the frequency response operator of system |7]) 



7» = u(iu£ - ry x g. 



(8) 



For a stable system Q, T{uj) describes the steady-state response to harmonic input signals across the 
temporal frequency uj. Namely, if the input is harmonic in i, i.e., 

d(y,t) = d^u,)^ 4 , 

with d(-,w) denoting a square-integrable spatial distribution in y, then the output is also harmonic in t with 
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the same frequency but with a modified amplitude and phase 



<p(y,t) = ([THd(, w )]b))e w = 0(y,Lj), 



The amplitude and phase of the output at the frequency oj are precisely determined by the frequency response 
operator T(u>), with 7ker denoting the kernel representation of the operator T. 

For the class of systems that we consider, the kernel representation of the frequency response operator 
7ker( • j • ; is a bounded matrix valued function on [a, b] x [a, b]. This implies that the operator T(lo) can 
be represented using the singular value (i.e., Schmidt) decomposition pT] . 

oo 

<f(y,u) = T(w)d(-,w) (y) = ^2 0-n(w)u„(y,w) (v n ,dj , (9) 

n = l 

where (•, •) is the standard Li [a, b] inner product, 

(vi,v 2 ) = / vl(y) v 2 (y)dy, 

J a 

and v^y) is the complex-conjugate-transpose of the vector vi(y). In ([9]), {u„} and {v„} denote the left and 
the right singular functions of the operator T associated with the singular value o~ n . These are obtained from 
the eigenvalue decomposition of the operators TT* and T* T, 

[T(w)7*(w)u*(-,«)](i0 - ^(wjflnd/.w), 
[T*(uj)T(uj)v n (-,Lj)}(y) = ctI(uj) v n (y,w), 

where T* is the adjoint of the operator T. The singular values are positive numbers arranged in descending 
order, 

CTl > cr 2 > ••■ > 0, 

and they are determined by the square root of the non-zero eigenvalues of TT* (or T*T). On the other 
hand, the singular functions {u n } and {v ra } form the orthonormal bases for the spaces of square integrable 
functions that the output (p and the input d belong to. 

From |9| we see that the action of the operator T(w) on d(y, ui) is determined by the linear combination 
of the left singular functions {u rl }. The product between the singular values, er n , and the inner product of 
the input d and the right singular function v„, |v„, d^, yields the corresponding weights. Thus, for d = v m , 
the output is in the direction of u m and its energy is determined by cr^ 



m ' 



d(y,w) = v m (y,w) ^> <p(y,w) = a m (u) u m (y , u) , 

implying that at any frequency w the largest singular value <ti(u>) quantifies the largest energy of the output 
for unit energy inputs. This largest energy can be achieved by selecting d(y,uj) = vi(y,uj), and the most 
energetic spatial output profile resulting from the action of T(io) is given by (p{y,uj) = <Ti(w) Ui(y, u>). 

In linear dynamical systems, spectral decomposition of the dynamical generators is typically used to 
identify instability. Appearance of the eigenvalues with positive real part implies exponential temporal growth 
of infinitesimal fluctuations and the associated eigenfunctions characterize spatial patterns of these growing 
modes. For systems with normal dynamical generators (i.e., operators that commute with their adjoints) 
the eigenfunctions are mutually orthogonal and the eigenvalues provide complete information about system's 
response. However, for systems with non-normal generators eigenvalues may give misleading information 
about system's responses. Even in the stable regime, non-normality can cause (i) substantial transient 
growth of fluctuations before their asymptotic decay; (ii) significant amplification of ambient disturbances; 
and (iii) substantial decrease of stability margins. We note that singular value decomposition of the frequency 
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Figure 2: Block diagrams of (a) the frequency response operator T: d H ► <p; and (b) the adjoint operator T*: f >-» g. 



response operator represents an effective tool for capturing these non-modal aspects of the system's response. 

In what follows, we describe the procedure for reformulating the frequency response operator ^ into cor- 
responding two point boundary value problems that are given by either an input-output differential equation 
or by a spatial state-space representation. These can be solved with superior accuracy using recently devel- 
oped computational tools [H]. We illustrate the utility of our developments on an example from viscoelastic 
fluid dynamics, where standard finite dimensional approximation techniques fail to produce reliable results. 

3. Two point boundary value representations of 7~, T* , and TT* 

In this section, we hrst describe the procedure for determining the two point boundary value representa- 
tions of the frequency response operator ^ . These are given by either a high-order input-output differential 
equation or by a system of first-order differential equations in spatial variable y. We then discuss the proce- 
dure for obtaining corresponding representations of the adjoint operator T* and the operator TT* ■ 



3.1. Representations of the frequency response operator T 

The application of the temporal Fourier transform to Q yields 



(\w£ - T) <p(y,u) 



H(j>(y,uj), 



(10a) 
(10b) 



where we have omitted hats from the Fourier transformed fields for notational convenience (a convention that 



we adopt from now on). System (10) represents an w-parameterized family of ordinary differential equations 
(ODEs) in y, with boundary conditions at a and b. From the definitions of the operators {£, Q, H.} 
described in Section 2.2 (10 1 can be represented by the following system of differential equations 



[AoH(y) = [B d](y), 
T: I <p(y) = [CofHv), 

= A/u0(t/), 



(11) 
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Here, DW<f>j = d l 4>j/dy l , E a and E& denote the point evaluation functionals at the boundaries, e.g., 

and {W a .i, Ws J are constant matrices that specify the boundary conditions on <p. For notational conve- 



nience we have omitted the dependence on ui in (11), which is a convention that we adopt from now on. 
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Here, n, m, k, and i denote the highest differential orders of the operators Aq, B , C , and A/"o, respectively. 
If the number of components in </>, d, and ip is given by s, r, and p, then {ai(y)} are matrices of size s x s 
with entries determined by the coefficients of the operator (iuj£ — J 7 ); {/3;(y)} are matrices of size s x r with 
entries determined by the coefficients of the operator Q\ and {~fi(y)} are matrices of size p x s with entries 
determined by the coefficients of the operator H. We also normalize the coefficient of the highest derivative 
of each fa to one, i.e., 

&7li,ii — 1; ^ — 1, . . . , 5, 

where a ni .a is the iith component of the matrix a ni , and rii identifies the highest derivative of fa. I n order 



to make sure that the input field d in ( 11 ) does not directly influence the boundary conditions and the output 



field (p, we impose the following technical assumptions on system (11), 



£ < n, m < n 



k < n 



This assumption is satisfied in most physical problems of interest. 

Alternatively we can bring ( |11[ ) into a system of first-order differential equations (in y). This can be done 
by introducing state variables, {xi(y)}, where each of the states represents a linear combination of <p and d, 
and their derivatives up to a certain order. A procedure for converting a high-order two point boundary value 



realization (11) with spatially varying coefficients to a system of first-order ODEs is described in Appendix 



[A") This transformation yields the spatial state-space representation of the frequency response operator T 



T : 




A (»x(y) 
Co(y)x(y), 

N a x(a) + 



+ B (y)d(y), 
N t x(6), 



(12) 



where x is the state vector, Ao, Bo, and Co are matrices with, in general, spatially varying entries, N a and 
N;, are constant matrices that specify the boundary conditions, and x' = dx/dy. To avoid redundancy in 
boundary conditions, N a and are chosen so that the matrix [ N a N& ] has a full row rank. We note 
that ( 12 ) is well-posed (that is, it has a unique solution for any input d) if and only if [35] 



dct(N a + N b * (M)) + 0, 



where ^o(y,r]) is the state transition matrix of Ao(y), 

d*o(y,?7) 



dy 



Ao(t/)* (y,»7), ^o(v,v) = I, 



and det (•) is the determinant of a given matrix. 

For the ID diffusion equation of Section [2.1.1| the input-output differential equation and the corresponding 
spatial state-space representation of the frequency response operator are given by (J3|and Q, respectively. 
Note that the boundary conditions (3b) can be rewritten into the form required by (11), 



3. 2. Representations of the adjoint operator T* 

We next describe the procedure for obtaining the two point boundary value representations of the adjoint 



of the frequency response operator, T, T*: f n- g; see figure 2(b) As shown above, the operator T can be 
recast into the input-output differential equation ( |11[ ), and the corresponding representation of T* is given 
by 

' [A&xl>](v) - [C* f](y), 
T : { g(y) = [BSVKiO, (13) 

= AftVO/). 
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Figure 3: A cascade connection of T* and T with TT*: 1 1-> tp. 

Here, the adjoint operators are [2Ql 123] 

n 

[^](y) - £(-ir[DW«v)](y), [Qf](y) = £ (-if [d« ( 7 |f)] 

m £ 

WKy) = £(-ir[D«(/3*V0l(y), = ^(W^E a + W^E b )[D«^(y), 



i = 



i = 



where a*, /3*, and 7* are the complex-conjugate-transposes of the matrices cc^, /3j, and The boundary 
conditions on the adjoint variable i/> are determined so that the boundary terms vanish when determining 
the adjoint of the operator Aq. A procedure describing how to determine the boundary conditions of the 
adjoint system is given in [2Ql Section 5.5]. 

On the other hand, the state-space representation of the adjoint of the operator T is given in [22] 



z'(y) = -A* (y)z(y) - C* (y)f(y), 
T*:l g(y) = B* (y)z(y), 

= M a z(a) + M b z(b), 



(14) 



where Aq, Bq, and Cq denote the complex-conjugate-transposes of the matrices A , B , and Co- The 
boundary condition matrices M a and M& are determined so that [ M a Mb ] has a full row rank and 



[ M a M b ] 



-N? 



= 0. 



(15) 



A procedure for selecting M a and M b that satisfy these two requirements is described in pH Section 3.1]. 
Furthermore, we note that the well-posedness of the adjoint representation ( fl4| is guaranteed by the well- 
posedness of T. 

For the ID diffusion equation of Section 2.1.1 the adjoint of the operator T(uj) described by ([3| has the 
following input-output representation 



(D^+itjl)^) = f(y), 
g(y) = -ip{y), 



(16) 
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As specified in (14 1, the state-space representation of T*{ui) is determined by taking the appropriate complex- 
conjugate-transposes of the corresponding matrices in Q with the following boundary condition matrices 



Mi 



3.3. Representations of TT* 

From the above described representations of T and T* , we can determine corresponding representations 
of the operator TT*: f i-> <p. As illustrated in figure [3j this operator represents a cascade connection of 
the frequency response operator T and its adjoint T* ■ The input-output differential equation for TT* is 
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obtained by equating the output of T* in (13) with the input of T in |TT] ), i.e., d = g, yielding 



[A£](y) = [Bf](y), 
TT : I <p{y) = [Ce](y), 
= Af$(y), 



(17) 



where 



<Kv) 



N = 



W * 



, A 



B = 



Ao 




-B Q B* 
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C=[C Q ] 



Similarly, the spatial state-space representation of TT* is obtained by equating the input d in (12 1 to the 



output g in (14), which yields 



q'(y) = A(y)q(y) + B(y)f(y), 
TT:\ <p(y) = C(y)d(y), 

= L a q(a) + L b q(6), 



(18) 



with 



q(y) 
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A (y) B (»)B5(y) 
-AS(y) 



C(y) = [ C (y) ] , 



L a = 

Since a cascade connection of two well-posed systems is well-posed, the existence and uniqueness of solu- 



N b 
M b 



tions of (17) and (18) is guaranteed by the well-posedness of the corresponding two point boundary value 
representations of T and T*. 

We next present a procedure for computing the largest singular value of T using the above representations 
of the operator TT*. 

4. Computation of the largest singular value of T 



In this section, we utilize the structure of the two point boundary value representations (17) and (|18[) of 



TT* to develop a method for computing the largest singular value of the frequency response operator T(w), 

^Lx(TM) = A max (TMT*M), 
where A max (-) denotes the largest eigenvalue of a given operator. In what follows, we present the procedure 



for computing the eigenvalues of TT* using both input-output (17) and state-space (18) representations 
of TT*. This is done by first recasting the system of differential equations into a corresponding integral 
formulation; we then employ recently developed automatic Chebyshev spectral collocation method [5] to 
solve the eigenvalue problem for the resulting integral equation. Note that the eigenfunction corresponding 
to the largest singular value identifies the output of the system that is most amplified in the presence of 
disturbances. Similar procedure can be used to determine the principal eigenfunction of the operator T*T, 
thereby yielding the input that has the largest influence on the system's output. 



The solution to a two point boundary value problem (17) can be obtained numerically by approximating 
the differential operators using, e.g., a pseudo-spectral collocation technique [TJ [2 [3j 0]. For differential 
equations of a high-order, the resulting finite-dimensional approximations may be poorly conditioned. This 
difficulty can be overcome by converting a high-order differential equation into a corresponding integral equa- 
tion |25j . This conversion utilizes indefinite integration operators that are characterized by condition numbers 
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that remain bounded upon discretization refinement. The procedure for achieving this conversion, described 
in Section 4.2 extends the result of (55] from a scalar case to a system of high-order differential equations. 
Furthermore, in Section 4.3 we show how a spatial state-space representation ( 18 1 can be transformed to an 
equivalent integral form. We then employ Chebfun's function eigs to perform the eigenvalue decomposition 
of the resulting system of equations. 

4-1. An illustrative example 

We first illustrate the procedure for converting a differential equation into its corresponding integral form 
using the ID diffusion equation ([3]), 



(d^ - iwl) <f>(y) = -d(y), 



Ei 4>{y) 



System ( 19 1 can be converted into an equivalent integral equation by introducing an auxiliary variable 



Integration of ( 20 ) yields 



<t>(y) 



'12 



(y)- 



(19a) 
(19b) 

(20) 



GO + h, 



v(r)i) dru dr) 2 + k x (y + 1) + k 2 



(21) 



[J< a > v] (y) + k , 



where jW and j' 2 ' denote the indefinite integration operators of degrees one and two, the vector k = 
[ k 2 ki ] contains the constants of integration which are to be determined from the boundary condi- 
tions ( 19b ), and 

= [ 1 (y + 1) ]. 



The integral form of the ID diffusion equation is obtained by substituting (21 1 into (191, 



(7 - iuj J {2) ^jv{y) - iwif (2) k = -d(y), 



1 
1 2 



A: 2 
ki 



Now, by observing that 



E-i 



(y) 



Ei 



(y) 



(22a) 
(22b) 



v(rf) drj = 0, 



we can use ( 22b ) to express the constants of integration k in terms of v 

(y) - 



k 2 ' 


1 


2 


" 




" " 


h _ 


~2 


-1 


1 




1 



J< 2 





-1/2 



Ei 



(!/)• 



(23) 



Finally, substitution of (23) into (22a) yields an equation for v, 

1 1 - iuJ™ + \w(y + 1) Et \v{y) 



<y)- 



(24) 



Invertibility of the matrix that multiplies the integration constants k = [ k 2 fci ] T in ( 22b ) facilitated 
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derivation of an explicit expression for k in terms of v. For the ID reaction-diffusion equation with homoge- 
nous Neumann boundary conditions, 



(D^ -cl - \u3l) cf>(y) = - d(y) 
£>«0l (y) = 



E, 



substitution of (21 1 to (25) yields 

(i - (iw + c) J (2) ) v{y) - (iw + c) k = - d(y), 



1 
1 



k 2 



(y) = 



(25a) 
(25b) 

(26a) 
(26b) 



A positive reaction rate c in (25a) ensures stability in the presence of Neumann boundary conditions. 



Lack of invertibility of the matrix that multiplies the integration constants in (26b I is an obstacle to 
determining k explicitly in terms of v. Instead, the dependence of v on k and d can be obtained from (26a 

(y) = (i - (iw + c) J^)' 1 ((iw + c)K^k - d(yj) 



Now, substitution of (27) to (26b) yields 

k = G 

where the matrix G is given by 
G 



E 1 J« (i ~ Qu + c) J^) d(y), 



(27) 



(28) 



" 1 " 




" " 


1 


+ 


1 



Ei J (1 > (i - (iw + c) j( 2 )) 1 (iw + c) A^. 



Finally, an equation for z/ is obtained by substituting (|28|) into (|26a 
(/ - (iw + c) J (2) ) = { (iw + c)K^G- 1 



ExW (i - Qw + c)J®) 1 - A d(y). (29) 



Systems ( 24 ) and ( 29 1 only contain indefinite integration operators and point evaluation functionals which 
are known to be well-conditioned. This is a major advantage compared to their corresponding input-output 
differential equations (19) and (25). 

4-2. Integral form of a system of high- order differential equations 



We now present the procedure for converting a system of high-order differential equations (17), 

' [A£](y) = [Bf](y), 
rr : I <?{y) = [C£](y), 
= Af$(y), 



(30) 



to an equivalent integral form. The input and output vectors f(j/) and <p(y) have p elements, £(?/) is a 
2s-vector, and the operators in (30) are given by 



A = ^ ai (t/)D« B = £bi(|/)DW, C = ^Cj^JDW, AA = ^](Y a , 4 E a + Y M E b )D 



(i) 



i = 



i = 



i = 



i = 



12 



As illustrated in Section 4.1 instead of trying to find the solution £ to |17| ) directly, we introduce two auxiliary 
variables, v and k. The ith component of the vector u(y) = [ v\(y) ... V2s{y) ] T is determined by 



(y), 



(31) 



where rii denotes the highest derivative of Zi in 

[At](y) = [Bf](y). 



Integration of (31) yields 



(y) = \j (ni - j) ^] (y) + K^-j) k,-, j = 0, 



(32) 



where £ C ni is the vector of integration constants which are to be determined from the boundary condi- 
tions, is the indefinite integration operator of degree with = 0, and K^ ni > is the matrix with 
columns that span the vector space of polynomials of degree less than m, 

KM = [K (y) K x {y) ■■■ K ni _ x {y) ], = 0, 

K (y) = 1, K,(y) = ^(y-a)\ j > 1. 



Substitution of (32) into (30) yields the integral representation of the operator TT* , 



TT : ^ 





' C n 


£12 ' 




V 




' B ' 




£21 


£-22 




k 








<P = [ Vi V2 



V 

k 



where 



C n = ^a,(y)j("- 4 ), £ X2 = ^ ai (y)K("- 1 ), 

i=0 i=0 

e e 

£21 = ^Y M E b j("- l \ C22 = £(Y 0| iE + Y M E 6 )K( n 

Pi 



i = 
k 



-i) 



i = 
k 



E c ^(y) J( "" l) ' 

i = 

j(ni-i) 



Vn 



j(n-i) _ 



2f(ni-») 



i = 



j{n 2 s-i) 

jW = 0, #W =0, i < 0. 
Using (33) we can determine an expression for the integration constants, 

£22 k = - [£ 2 i 1/] (y). 



If the matrix £22 is invertible, equation (34) in conjunction with (33) yields 

"(V) = [(Ai - Aa^^ai)" 1 (Bf)] (»), 
= [(^1 - Pa £22* £21 H (y), 



(33) 



(34) 

(35a) 
(35b) 
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and the representation of the operator TT* is obtained by substituting (35a) into (35b I. Thus, determination 
of the left singular functions {u„} of the operator T amounts to solving the following eigenvalue problem 



(Vi - V 2 £^£ 21 ) (Cn - £ 12 £^£ 21 ) (Bu n ) (y) = a 2 n u n (y), 



(36) 



where o~ n denotes the corresponding singular value of T . 



On the other hand, if £22 is singular, we can determine an expression for v in terms of k and f from (33), 

u(y) = [C^Bi] (y) - Crf £ 12 k. (37) 



Furthermore, substitution of ( |37| into (|34j) yields 

k = -G _1 £2i [C^Bf] (y), 

where the matrix G is given by 

G = £ 22 — £ 2 i £u £y 2 . 
This expression for k in conjunction with ( |33[ ) yields 

v{y) = [Crf (B + £ 12 G 1 £ 21 £^B) f] (y), 
<f(v) = [Piv](y) - [V 2 G-' £ 21 £^Bi\ (y). 



(38) 



(39a) 
(39b) 



The integral representation of the operator TT* can be obtained by substituting (39a) into (39b), and the left 
singular pair (o~ n ,u n ) of the operator T is determined from the solution to the following eigenvalue problem 



[(7>i C£ + V x £~il £ 12 G 1 £21 £^1 - V 2 G 1 £ 21 C£) (8u„)] (y) = a\ u n (y). 
4-3. Integral form of a spatial state-space representation 



(40) 



We next describe a procedure for transforming a spatial state-space representation (18), 

q'(y) = A(y)q(y) + B(y)f(y), 
TT- : { <p(y) = C(y) q (y), 

= L a q(a) + L b q(6), 



(41) 



into a system of first-order integral equations. In a similar manner as in Section |4.2[ we introduce two 
auxiliary variables u and k so that 



»{y) = q'(y) q(y) = [Jv] (») + k, 

where J is a block diagonal matrix of the first order indefinite integration operators 



(42) 



Substitution of (42) into (41 ) yields a system of first order integral equations for the operator TT*, 

u{y) = A(y) [Jv](y) + A(y) k + B(y) f(y), (43a) 
<p(y) = C(y) [J u] (y) + C(y)k, (43b) 
0= (L a E Q + L 6 E 6 )[Ji/](y) + (L„ + L 6 )k. (43c) 
An expression for v in terms of the forcing f and the integration constants k can be obtained from (43a), 



(I-AJr^Bf) (y)+ (I -A J)" 1 A (y)k. 



(44) 
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Figure 4: Channel flow geometry. 



Furthermore, substitution of (44 1 into (43c) yields 



k = H 1 L b E fc 



J(I-AJ)- 1 Bf (y) 



(45) 



where H is a matrix given by 



H 



Lb Eb 



3 (I- Ajy 1 A] (y) + L a + L 



Finally, substitution of (|44|) and (|45|) into (|43b|) yields 



CJ(I-AJ)~ 1 Bf (y) - CH ^LtEhJ^ AJ) _1 Bf 



(?y) 



(46) 



- CJ(I AJ) 1 AH ^LbEbJ^ AJ) x Bf (y), 
where invertibility of the matrix H follows from the well-posedness of the two-point boundary value prob 



lem (41 1. Thus, the singular values a n and the associated left singular functions u„ of T can be obtained by 



solving the following eigenvalue problem 



CJ(I - AJl^Bu, 



(y) 



CH^LiEiJfl -AJ) _1 Bu r , 



(y) 



CJ(I A3) 1 AH 1 L b E b 3 (I AJ) _1 Bu 



(v) = ^u„(y). 



(47) 



In summary, the principal left singular pair of the operator T can be determined by rewriting either the 



input-output differential equation (17) or the system of first-order differential equations (18) representing 



TT* into their respective integral forms ( 33 1 and ( 43 ) . The resulting eigenvalue problems ( 36 ) and ( 47 1 



are solved using Chebfun [pj. The detailed discussion on how Chebfun can be used to solve the eigenvalue 



problems ( 36 ) and ( 47 ) is relegated to Appendix B 



5. Examples 

We next use our method to study frequency responses of two systems from fluid mechanics: the three- 
dimensional incompressible channel flow of Newtonian fluids, and the two-dimensional inertialess channel 
flow of viscoelastic fluids. In the latter example, we show how numerical instabilities encountered when using 
finite dimensional approximation techniques can be alleviated. The utility of theoretical and computational 
tools of this paper goes beyond fluids; they can be used to examine dynamics of a broad class of physical 
systems with normal or non-normal dynamical generators, and spatially constant or varying coefficients. 



5.1. Three-dimensional incompressible channel flows of Newtonian fluids 

We first study the dynamics of infinitesimal three-dimensional fluctuations in a pressure-driven channel 
flow with base velocity U(y) = l — y 2 ; see figure [4] for geometry. As shown in [27], the linearized Navier-Stokes 
(NS) equations can be brought to the evolution form |t| with state = [ <f>\ <p2 ] T , where 4>i an d 4>2 are 

the normal velocity and vorticity fluctuations. Furthermore, d = [ d\ d2 c?3 ] T and cp = [ u v w ] T 
are the input and output fields whose components represent the body forcing and velocity fluctuations in the 
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CO 

Figure 5: Two largest singular values of the frequency response operator for the linearized Navier-Stokes equations as a function 
of the temporal frequency ui in the flow with R = 2000, k x = 1, and k z = 1: blue X, <ti(T); and red o, o"2(T). 





(b) 

Figure 6: Streamwise velocity fluctuation development for largest singular value of the frequency response operator in a pressure- 
driven channel flow with R = 2000, k x = k z = 1, (a) uj = —0.385, and (b) ui = —0.982. High and low velocity regions are 
represented by red and green colors. Isosurfaces of u are taken at ±0.55. 



three spatial directions, x, y, and z. Owing to translational invariance in x and z, Q is parameterized by the 
corresponding wave numbers k x and k z with the boundary conditions on the normal velocity and vorticity, 



Mkx,±l,k z ,t) = DWMk x ,±l,k z ,t) = 0, 

4, 2 (k x ,±l,k z ,t) = 0, k x ,k z eR, t>0. 



The operators in Q are given in Appendix C and, for any pair of k x and k z , they are matrices of differential 
operators in y G [—1, 1]. 

In what follows, we set the Reynolds number to R = 2000, k x = k z = 1 and compute the singular values 
of T using the method developed in Section 4.2 Figure [5] shows two largest singular values, a± and (72, of 
the frequency response operator T for the linearized NS equations as a function of the temporal frequency 
uj. The largest singular value a\ exhibits two distinct peaks at w ~ —1 and uj ~ —0.4. These peaks are 
caused by different physical mechanisms which can be uncovered by investigating responses from individual 
forcing to individual velocity components |27j . The discussion of these mechanisms is beyond the scope of 
this paper. 

Figure [6] shows the isosurface plots of the most amplified streamwise velocity fluctuations corresponding 
to the two peaks shown in figure[5j These structures are purely harmonic in x, z, and t, and their profiles in y 
are determined by the the left principal singular functions of the frequency response operator at uj = —0.385 
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z z 

(a) (b) 



Figure 7: Streamwise velocity (color plots) and vorticity, u> x = d y w — d z v, (contour lines) fluctuation development for largest 
singular value of the frequency response operator in the cross section of a pressure-driven channel flow with R = 2000, k x = 
k z = 1, (a) ui = —0.385, and (b) uj = —0.982. Red color represents high speed and blue color represents low speed streaks. 



and uj = —0.982. For uj = —0.385, u is localized in the near-wall region. On the other hand, for uj = —0.982 
the fluctuations occupy the center of the channel. The development of the streamwise velocity (color plots), 
and streamwise vorticity uj x = d y w — d z v (contour lines) fluctuations in the channel's cross-section is shown in 
figure]?] For uj — —0.385, the most amplified set of fluctuations results in pairs of counter rotating streamwise 
vortices that generate high and low velocity in the vicinity of the lower and upper walls. In contrast, for 
uj = —0.982 there is a large concentration of arrays of counter rotating streamwise vortices in the center of 
the channel. Even though the spatial patterns identified by our analysis represent an idealized view of the 
flow, their utility in understanding the early stages of transition to turbulence has been well-documented [11] . 



5. 2. Inertialess channel flow of viscoelastic fluids 

We next compute the frequency responses of the inertialess flow of viscoelastic fluids presented in Sec- 
tion [2?L2] This example illustrates the utility of our method in situations where standard finite dimensional 
approximations may fail to produce accurate results. For this example, the input-output and spatial state- 
space representations of the frequency response operator are given in Appendix D We compute the largest 
singular value using the procedure described in Section [4] and provide comparison of our results with those 
obtained using a pseudo-spectral collocation method [18] . 

It is well-known that inertialess flows of viscoelastic fluids exhibit spurious numerical instabilities at high- 
Weisssenber numbers [TBI [17]. I n view of this, we fix k x — 1, = 0.5, and uj = and examine the effects of the 
Weissenberg number, We, on the frequency response. We first compute the largest singular value of T using 
a pseudo-spectral collocation method [TH]. This is achieved by approximating the operators in the input- 
output representation (17) of TT* with differentiation matrices of different sizes. Figure 8(a) shows that 
°max converges as the number of collocation points, N, increases from 50 to 200 for 1 < We < 9. However, 
for We > 9 the increased number of collocation points in y does not necessarily produce convergent results; 
see figure [8(b) [ Furthermore, in certain cases, the eigenvalues of the operator TT* computed using pseudo- 
spectral method have large negative values. This is clearly at odds with the fact that TT* is a non-negative 
self-adjoint operator, which indicates that the negative eigenvalues arise from numerical artifacts. 

Figure s | 8(c)| and |8(d)| show the largest singular value of the operator T computed using the method 
of Section |4[ For 1 < We < 9, the largest singular values obtained in Chebfun for both input-output and 
spatial state-space integral representations of TT* are equal to each other and they agree with the results of 
pseudo-spectral method; see figure 8(c) For We > 9 we see that the largest singular value computed using 



Chebfun exhibits nice trends as We increases. Furthermore, automatic Chebyshev spectral collocation method 
employed by Chebfun makes sure that grid point convergence of the singular values is satisfied. We note that 
the singular values computed using the input-output and spatial state-space integral representations of TT* 
are equal to each other for We < 12. On the laptop used for computations, Matlab has experienced memory 
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Figure 8: The largest singular values of the frequency response operator for an inertialess shear-driven channel flow of viscoelastic 
fluids as a function of We at k x = 1, f3 = 0.5, and ui = 0. Results are obtained using: (a) and (b) Pseudo-spectral method with 
N = 100, blue o; N = 150, red *; and N = 200, green □; (c) and (d) Chebfun with integral forms of input-output differential 
equations, blue A; and spatial state-space representations, red V. 
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Figure 9: Principal singular functions of the frequency response operator for inertialess shear-driven flow of viscoelastic fluids 
with We = 20, k x = 1, and = 0.5. First column: real part of « m ax; second column: imaginary part of u max . Results are 
obtained using: (a) and (b) Pseudo-spectral method with N = 50, red X; N = 100, blue o; TV = 200, green □; (c) and (d) 
Chebfun with integral form of input-output differential equations. 



issues when solving the eigenvalue problem in the state-space formulation (47) for We > 12. These memory 
issues may arise from solving a large system of linear equations internally in Chebfun. We have not observed 
any problems with memory when solving eigenvalue problem obtained from the the input-output equation 
in its integral form ( 36 ) . We further note that the singular values can be computed accurately using the 
input-output integral representation at much higher Weissenberg numbers. 

We next present the principal singular functions corresponding to the streamwise and normal velocity 
fluctuations in a flow with We — 20. These are obtained using pseudo-spectral method and Chebfun with the 
input-output integral representation. Figures |9(a)| and |9(b)| show the spatial profiles of velocity fluctuations 
that experience the largest amplification in the presence of disturbances. These profiles are obtained using 
pseudo-spectral method with different number of collocation points. Note the lack of convergence as the 
number of collocation points is increased. On the other hand, Chebfun does not suffer from numerical 
instabilities, and the corresponding principal singular functions exhibit nice symmetry with respect to the 
center of the channel; see figures 9(c) and |9(d)| Similar trends are observed for larger values of We. 



19 



6. Concluding remarks 



We have developed a method for computing the principal singular value and the corresponding singular 
functions of the frequency response operator for distributed systems with a spatial variable that belongs to 
a compact interval. Our method avoids the need for numerical approximation of differential operators in the 
evolution equation. This is achieved by recasting the frequency response operator as a two point boundary 
value problem, a formulation well-suited for employing Chebfun computing environment. When dealing with 
spatial differential operators of high order our method exhibits two advantages over conventional techniques: 
numerical ill-conditioning associated with high-order differential matrices is overcome; and boundary condi- 
tions are easily implemented and satisfied. We have provided examples from Newtonian and viscoelastic fluid 
dynamics to illustrate the utility of our developments. 

Our method has been enhanced by the development of easy-to-use Matlab functions which take the 
system's coefficients and boundary condition matrices as inputs and yield the desired number of left (or 
right) singular pairs as the output. The coefficients and boundary conditions of the adjoint systems are 
automatically implemented within the code using the method described in this paper. The burden of finding 
the adjoint operators and boundary conditions is thus removed from the user who can instead focus on 
interpreting results and understanding the essential physics. 

Even though we have confined our attention to computation of the frequency responses for PDEs, the 
developed framework allows users to employ Chebfun as a tool for determining singular value decomposition 
of compact operators that admit two point boundary value representations. In particular, our approach paves 
the way for overloading Matlab's command svds, from matrices to compact operators. 

While the body of the paper focuses on PDEs with distributed input and output fields, by considering 
an Euler-Bernoulli beam with boundary actuation in |Appendix E[ we illustrate how Chebfun can be used to 
compute frequency responses of systems with boundary inputs. This problem turns out to be much simpler 
than the problems with distributed inputs, and it can be implemented with only few lines of code in Chebfun. 
We also use this example to demonstrate the utility of integral formulation in producing accurate results 
even for systems with poorly scaled coefficients. 

In all examples that we considered, it is much more efficient to compute the eigenvalue pairs for a system 



of high-order integral equations (36) than for a system of first-order integral equations (47). We believe that 
larger number of dependent variables is reducing efficiency of computations that rely on spatial state-space 
representation. We note that Chebfun automatically adjust the number of collocation points in order to 
obtain solutions with an a priori specified tolerance. The computational speed can be increased by lowering 
this tolerance using the following command in Matlab 

chebf unpref ( ' res ' , tolerance) . 

Our ongoing efforts are focused on employing Chebfun as a tool for computing the peak (over temporal 
frequency) of the largest singular value of the frequency response operator. In systems and controls literature, 
sup^ <7 max (7~(w)) is known as the Hoo norm and its computation requires identification of purely imaginary 
eigenvalues of a Hamiltonian operator in conjunction with bisection [28] . In addition to quantifying the worst- 
case amplification of purely harmonic (in time) deterministic (in space) disturbances, the inverse of the 
norm determines the size of an unstructured modeling uncertainty that can destabilize the nominal system. 
Thus, large frequency response peaks indicate small stability margins (i.e., poor robustness properties to 
modeling imperfections) , and they are a reliable predictor of systems in which small modeling imperfections 
can introduce instability. This interpretation of the norm is closely related to the notion of pseudospectra 
of linear operators [7] and it has been used to provide useful insight into dynamics of systems with non-normal 
generators [51 HT]. 
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Supplementary material 

All Matlab codes for computing frequency responses are available for download at 

www.umn.edu/^mihailo/software/chebfun-svd/ 

Appendix A. Conversion to a spatial state-space realization 

We next describe how a high-order ODE with spatially varying coefficients can be converted to a family 



of first-order ODEs (12 1. We consider the following ordinary differential equation with boundary conditions 

n — 1 rci 

0(»)(y) = 



- OiM^iy) + Y,0i(y)dW(y), m <n-£, 

i = Q 

k 

v{v) = ^2in{v)4>®{v), k < 

i = 

e 

= £jVi, o 0W(a) + N itb <t>®(b), I < n, 



» = 



n — to, 



(A.la) 
(A.lb) 
(A.lc) 



i = Q 



where <jr- % > — d 1 <fi / 'dy l . Since coefficients in (A.la) are spatially varying, the standard observer and 

controller canonical forms cannot be used to obtain a system of first-order ODEs (12 1. Instead, we introduce 
a new variable w(y), 



i = 



and substitute (A.2) into (A.la) to obtain 

n- 1 



» = 



Here, a state-space realization of (A.3) is given by the controller canonical form, 

z'(y) = Ai(j/)z(y) + e n w(y), 

<f>(y) = efz(y), 

where 

Ax(y) = 



(A.2) 



(A.3) 



(A.4a) 
(A.4b) 



(n-l)xl 


I (n-l)x(n-l) 


-u (y) 


-ai(y) ••• -a„_ 1 (y) 



and e; is the ith unit vector. It is a standard fact that the solution to (A.4) is given by 



z(?y) = *l(y,a)z(a) + / »/)e„w(77)dTO 

*/ a 



where &i(y,r)) is the state-transition matrix of Ai(y). Substituting (A.2) into (A. 5) yields 
z(y) = *i(y,a)z(a) 



(A.5) 



(A.6) 
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Application of integration by parts to the integral in (A.6) along with a change of variables leads to the 
following two point boundary value state-space representation of ( A.l[ ) 



where 



x (y) 

A (y) 
C (tf) 



x'(y) = A (y)x(y) + B d(y), 
<P(V) = C x(y), 

= N a x(a) + N b x(6), 



m — 1 I rn — i 

<y) -EE Qi-i(A-+i(v)) ] <*<%), 

i=0 \J=1 

m 

A x (y), B (tf) = ^Qi(ft(y)), 
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(A.7b) 
(A.7c) 
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We note that, for a given function /3, Qi can be recursively determined from 

d 



QMy)) = AtMCU-tMv)) 
Qo(P(y)) = e„ /%). 



dy 



Qi-i(/3(2/)), i = 1, m, 



Appendix B. Implementation of eigenvalue problems in integral formulation using Chebfun 



The eigenvalue problems (36) and (47) derived in Sections 4.2 and 4.3 are solved using Chebfun. Here, 



we show how to implement the functions and operators in Chebfun to solve (36); a similar procedure can 



be used to solve (|47j) . The eigenvalue problem ( 36 1 requires the construction of a number of operators and 



quasimatrices (terminology used by the authors of Chebfun to denote vectors of functions). The operator A 



in (30) is represented by the coefficients a;(y) which are functions determining columns of a quasimatrix. For 



example, consider the differential equations representing the operator TT* for the ID diffusion equation 

/(«), 



- yjjl -I 




' " 




' o " 


+ vjjI 








i 



m = [ i o ] 



' 1 " 




' Ci(-i) " 




" " 








" o " 







. _ 


+ 


1 











" 1 " 




' 6(-i) " 




" " 








" " 









+ 


1 












(B.l) 



The code used to generate operator A for the ID diffusion equation is given by 

%7o Operator A for the ID diffusion equation 
dom = domain(-l , 1) ; % domain of functions 

fone = chebf un(l ,dom) ; % fone(y) = 1 
fzero = chebf un(0, dom) ; % fzero(y) = 



22 



7, w is the temporal frequency and li is the imaginary unit 
7, (1,1) element of operator A 

All = [-li*w*fone, fzero, f one] ; % -i*w*xi_l + 0*D~{(l)}*xi_l + l*D~{(2)}*xi_l 
7, (1,2) element of operator A 

A12 = [-fone, fzero, fzero]; 7. -l*xi_2 + 0*D~{(l)}*xi_2 + 0*D~{(2)}*xi_2 
% (2,1) element of operator A 

A21 = [fzero, fzero, fzero]; 7. 0*xi_l + 0*D~{ (1) }*xi_l + 0*D~{(2)}*xi_l 
7 (2,2) element of operator A 

A22 = [li*w*fone, fzero, fone]; 7. i*w*xi_l + 0*D~{ (1) }*xi_l + l*D~{(2)}*xi_l 

7» form operator A using cell-array construction 
A = {All, A12; A21, A22}; 

The variable dom denotes the domain of the functions, and fone and fzero represent unit and zero functions. 
The dimension of each Chebfun's function in Matlab is oo x 1, where the first index represents the continuous 
variable y. Hence, the quasimatrices All, A12, A21, and A22 have dimensions oo x 3. Since the dimension 
of quasimatrices prohibits the construction of matrix of functions, we instead utilize Matlab's cell arrays 
(using curly brackets) to represent the operator A. The boundary condition matrices are given by 

Yal = [1, 0; 0, 0] ; Ya2 = [1 , ; , 0] ; 
Ybl = [0, 0; 1, 0]; Yb2 = [0, 0; 1, 0]; 
Ya = {Yal; Ya2>; Yb = {Ybl; Yb2> ; 

The code used to generate the quasimatrix K'™) is given by 

n = size(A,l); % number of states in your system of ODEs 

% determine the highest differential order of each component of \xi in the equations 
ni = zeros(n, 1) ; 
for j = l:n 

ni(j) = size( A{j,j>, 2) - 1; 

end 

% indefinite integration operator 
J = cumsum(dd) ; 

'/,% Construct each component of K 
Ki = chebfun(l ,dd) ; 
for j =2 : max(ni) 

Ki(: ,j) = J*Ki(: , j-1) ; 

end 

7. construct quasimatrix K using cell-array 
for j = l:n 

K{j} = Ki(: , l:ni(j)) ; 

end 

The indefinite integration operator is obtained using Chebfun's command cumsum. The variable ni contains 
the highest differential order of each state £j in the system. We next determine the matrix £22 appearing 
in (33 1 by applying the boundary condition operator Af to K. The following code is used to generate £22 

'/'/, Determine the matrix L_{22} 

7o loop through each component of \xi 

for j = l:n 

7« quasimatrix K associated with \xi_{j} 

Kj = K{j>; 

L22{j> = Ya{j} + Yb{j}*toeplitz([l zeros(l, ni(j)-l)], Kj ( b, : )); 

end 

The qausimatrix £12 is obtained by multiplying coefficients of the operator A with the quasimatrix K, 



23 



%7o Determine the functional operator L_{12} 

°/ loop through each component of L_{12}, which has size n x n 
for i = l:n 
for j = l:n 

'/, initialize the (i,j) component of L_{12} and 
7. get the quasimatrix K associated with \xi_{j} 
L12ij = 0; Kj = K{j}; 

7. get the (i,j) component of operator A 

Aij = A{i,j}; 

for indni = 1 : ni(j) 

L12ij = L12ij + diag( Ai j ( : , ind) )*Kj; 

Kj = [ chebfun(0,dd) , K j ( : , l:ni(j) - 1) ] ; 

end 

L12{i, j} = L12ij; 

end 

end 



The operator L\\ in (33) is realized using the following Matlab's commands 
%7o Determine the operator L_{11} 

7o loop through each component of L_{11}, which has size n x n 
for i = 1 : n 
for j = 1 : n 
7o get the (i,j) component of A 
Aij = A{i,j>; 

7o initialize (i,j) component of Lll with Aij_0 
Lllij = diag( Aij(:,l) ); 
for indni = 1 : ni(j) - 1 

Lllij = Lllij*J + diag( Ai j ( : , indni + 1) ) ; 

end 

Lllij = Lllij*J + diag( Ai j ( : , ni(j) + 1) ) ; 
Lll{i,j> = Lllij; 

end 

end 

The boundary point evaluation functional Ef, is easily constructed by 
Eb = linop(@(n) [zeros (1 ,n-l) 1], @(u) feval(u,b), dd) ; 

In a similar manner, the operator £21 is realized by 
'/'/, Determine the operator L_{21} 

7o loop through each component of L_{21} which has size of n x 1 
for j = l:n 

7 get the j component of the boundary condition matrix Yb 

Ybj = Yb{j>; 

L21j = Ybj (: ,l)*Eb; 

for indni =1 : ni(j) - 1 

L21j = L21j*J + Ybj(:, ind+l)*Eb; 

end 

L21{j> = L21j*J; 

end 



We note that the operators V\ and Vi in (33) can be constructed using similar procedure. We have shown 
how to construct all operators and quasimatrices appearing in (33). However, the eigenvalue problem (36) 
requires the operator £ 12 £22 £21 • This operator can only be realized using explicit construction [2 6) because 
Chebfun syntax does not allow this expression to be formed directly. 
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%% determining the operator H = L_{12> L_{22}~{-1> L_{21> 
°/ looping through each component of H which has size of n x n 
for i = l:n 
for j = l:n 

L12ij = L12{i,j>; 

L22j = L22{j>; 

L21j = L21{j>; 

°/o m-by-m discretization of H (discretized form) 

mat = @(m) L12ij ( chebpts(m,dom) , : )*( L22j \ L21j (m) ); 

% functional expression of H (functional form) 

op = @(v) L12ij*( L22j \ (L21j*v) ); 

7. explicit construction of a linear operator in Chebfun 
H{i,j} = linop(mat,op,dom) ; 

end 



end 

A similar procedure is used to construct the operator Vi £22 £21 ■ Finally, Chebfun's eigenvalue solver 
(eigs) is used to compute the eigenvalues and eigenfunctions. We note that we use similar method to construct 
the operators for the spatial state-space representation of the eigenvalue problem discussed in Section |4~3) For 
brevity, they are not presented here. All codes for solving the eigenvalue problems in the integral formulation 
using Chebfun are available at www.umn.edu/~mihailo/software/chebfun-svd/. 



Appendix C. Representations of the frequency response operator for the linearized Navier- 
Stokes equations 

In this section, we provide the input-output and spatial state-space representations of the frequency 
response operator for the linearized NS equations. The input-output differential equations for the three- 
dimensional incompressible channel flow are given by 

f (a 4 D( 4 ) + a 2 (y)D( 2 ) + a (y)) <f>(y) = (biD« + b ) d(y), 

= ( Cl D« + c ) <My), (CI) 
k 0=((W_i,iE_ 1 +W w Ei)DW + (W_i,„E_i + Wi, o Ei))0( y ), 



T : <^ 



u 
v 
w 
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where 



a<t(y) = 
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The spatial state-space representation of T is obtained by rewriting (C.l I into a system of first-order differ- 
ential equations given by ( 12 ) with the following matrices 
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The input-output and state-space representations of the adjoint of the operator T can be determined using 
the procedure presented in Section [X2] 



Appendix D. Representations of the frequency response operator for the inertialess channel 
flow of viscoelastic fluids 

We next show how to formulate the input-output and spatial state-space representations of the frequency 
response operator for the inertialess flow of viscoelastic fluids. We begin by rewriting pi into the input-output 
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representation ( |11[ ), 

f + a 3 (y)£ (3) + « 2 (?y)^ (2) + ai^)^ 1 ) + ao(v))V(l/) = (My)D (1) + b (y)) d(y), 



T : < 

where 
a o(y) 

«i(y) 
«2(y) 
«3(y) 
My) 

Cl 



= (ci£>W + co)V(y), 
. 0= ((W-lxB-i + W M £?i)i?« + CW_i, S_i + W l!0 ia))V(y), 
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The spatial state-space representation of T is obtained by rewriting (D.l ) into a system of first-order differ 
ential equations. Using the procedure described in | Appendix A| yields 
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The input-output and state-space representations of the adjoint of the operator T can be determined using 
the procedure described in Section |X2) 



Appendix E. Frequency response of an Euler-Bernoulli beam 

In this section, we consider an Euler-Bernoulli beam that is clamped at the left end and subject to a 



boundary actuation u(t) at the other end; see figure E.10 for an illustration. The vertical displacement of 
the beam 4>(y,t) is governed by [2"§] . 



a Ei Ei 

^4>tt{V,t) + -^<t>tyy V y(y,t) + -J^<Pyyyy[V,t) = 0, y£[0,l], 



a Ei 

<t>yy{l,t) = 0, —^-(j)tyyy{l,t) 



0(O,i) = <f>y(Q,t) = 0, 

Ei 



(E.la) 
(E.lb) 

(E.lc) 



Here, the input u(t) denotes the force acting on the tip of the beam, I is the length of the beam, [i is the 
mass per unit length of the beam, Ei is the flexural stiffness, and a is the Voigt damping factor. 
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clamped 





forced by u(t 



Figure E.10: An Euler-Bernoulli beam that is clamped at the left end and subject to a boundary actuation at the other end. 



Equation (E.l) can be used to model the movement of a micro-cantilever in atomic force microscopy 



applications [30] with 



£ 
Ei 



240 x 10~ 6 m, n 
7.55 x 10~ 12 Nm 2 , a 



1.88 x 10" 7 kg/m, 
5 x 10~ 8 s. 



(E.2) 



In contrast to the body of the paper, the forcing u(t) does not enter to the equation as an additive input but 
as a boundary condition. We next show how easily frequency response in this case can be computed using 
Chebfun. 



Application of the temporal Fourier transform to (E.l I yields 

r Ei 



7» 



I 4 
0(0, u 



(1 + iua) <V'"{y,u) - mw 2 ^(2/,w) = 0, 



= 



<f>'(0,u) = 0, 
Ei 



(E.3) 



(1 + iua) <P"'{1,uj) = u(lo). 



At each uj, the mapping from u(ui) to (f)(y,uj) can be obtained by computing the solution to (E.3) with 
u(oj) = 1 using Chebfun. The energy of the beam is determined by 

E{lo) = i((0"(., W ),0"(.,c)> + c 2 (0(., W ),0(-, W ))), 

and it can be simply computed with the aid of Chebun's functions diff and cumsum. On the other hand, if the 
output is given by the vertical displacement at the right end of the beam, the frequency response is simply 



determined by the magnitude and phase of the complex number (/)(!, uj); see figure E.ll 



For parameters given by (E.2), even the use of Chebfun's differential operators to construct 



Ei 



+ icua)D {4) - Liu 2 I, 



with appropriate boundary conditions may lead to unfavorable conditioning of differentiation matrices. This 



can be alleviated by determining and solving instead the integral form of (E.3). The procedure for achieving 



this closely follows the method presented in Section |4~2) The Matlab code used for computing the frequency 
response with integral formulation can be found at www.umn.edu/^mihailo/software/chebfun-svd/. 
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